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Abstract 

The magnetic flux periodicity in superconducting loops is reviewed. Whereas quantization of 
the magnetic flux with hc/2e prevails in sufficiently thick loops with current free interior, the su- 
percurrent in narrow loops is either hc/2e or hc/e periodic with the external magnetic flux. The 
periodicity depends on the properties of the condensate state, in particular on the Doppler shift 
of the energy spectrum. For an s-wave superconductor in a loop with diameter larger than the 
coherence length the Doppler shift is small with respect to the energy gap, and the hc/2e 
periodic behavior of its flux dependent thermodynamic properties is maintained. However, for 
smaller s-wave loops and, more prominently, narrow (i-wave loops of any diameter R, the Doppler 
shift has a strong effect on the supercurrent carrying state; as a consequence, the fundamental 
flux periodicity is in fact hc/e. It is shown analytically and numerically that the hc/e periodic 
component in the supercurrent decays only algebraically as 1/i? for large d-wave loops. For nodal 
superconductors the discrete nature of the eigenergies close to the Fermi energy has to be respected 
in the evaluation of the Doppler shift. Furthermore, we investigate, whether the Doppler shift mod- 
ifies the supercurrent through Josephson junctions with (i-wave superconductors. For transparent 
junctions, the Josephson current behaves similar to the persistent supercurrent in a loop. These 
distinct physical phenomena can be compared, if the magnetic flux <!> = (/>• hc/e is identified with 
the phase variation of the order parameter Jy? through 27r(/> = 5ip/2. Correspondingly, the Joseph- 
son current can display a 47r periodicity in 5ip, if the Doppler shift is sufficiently strong which is 
true for transparent junctions of d-wave superconductors. Moreover, a An periodicity is also valid 
for the current-flux relation of field-threaded junctions. In the tunneling regime the microscopic 
theory reproduces the results of the Ginzburg-Landau description for sufficiently wide Josephson 
junctions. 
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I. INTRODUCTION 



The quantum mechanical wave function ip of particles moving in a multiply connected 
geometry has to be a unique function of the spatial coordinate. This condition leads to a 
discrete energy spectrum, because the phase difference of the wave function accumulated 
on a closed path has to be 2'iTk, where the integer k serves as a quantum number of the 
wave function. For a circular geometry, this phase winding number k represents the angular 
momentum hk of the particles. 

In the presence of a magnetic field B(r) = V x A(r), an additional term adds to the 
phase of the wave function: tp' = ip exp(— i 27r(e//ic) JJ"^ dr' ■ A(r')), where A(r) is the vector 
potential, e the charge of the electron, c the velocity of light, h is Planck's constant, and 
To an arbitrary space point within the system. The gauge transformed wave function if^' 
satisfies the Schrodinger equation with the vector potential A eliminated from the kinetic 
energy term. The new condition is that ip' acquires the phase factor exp(— z27r(e//ic)$) for a 
path C enclosing the magnetic flux $ = dr ■ A(r). This leads to a total phase difference of 
2Ti{k — e^/hc) on the closed path C. Because physical quantities are obtained by a thermal 
average over all possible k, they are periodic in $ with the fundamental period 

$0 = hc/e, (1) 

which is the flux quantum in the normal state. In particular, the persistent current </($) 
induced by the magnetic flux vanishes whenever $/<I>o is an integer. 

The effect described above is present in any system with sufficient phase coherence, and 
best known from the periodic resistance modulations of a microscopic metallic loop, pre- 
dicted first by Ehrenberg and Siday in 1948^ and in 1959 by Aharonov and Bohni^. Already 
ten years earlier, London predicted the manifestation of a similar effect in superconducting 
loops, where the phase coherence is naturally macroscopic^: the magnetic flux threading 
the loop is quantized in multiples of $o; because the interior of a superconductor has to be 
current free. London did not know about the existence of $o/2 flux quanta in superconduc- 
tors, but he already speculated that the supercurrent might be carried by pairs of electrons 
with charge 2e and that the superconducting flux quantum and hence the flux periodicity 
of the supercurrent is rather $o/2. This point of view became generally accepted after the 
'Theory of Superconductivity' by Bardeen, Cooper, and Schrieffer (BCS) was published in 
1957-. Direct measurements of magnetic flux quanta $o/2 trapped in superconducting rings 
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FIG. 1: Scheme of the pairing of angular-momentum eigenstates in a one dimensional metal loop 
for (a) = and (b) $ = $o/2, as used by Schrieffer ii^ to illustrate the origin of the $o/2 
periodicity in superconductors. Paired are states with equal energy, which leads to pairs with a 
center-of-mass angular momentum g = in (a) and g = 1 in (b) in units of h. 

followed in 1961 by Doll and Nabauer- and by Deaver and Fairbank-, corroborated later by 
the detection of $o/2 flux lines in the vortex phase of type II superconductors^*^. 

For thin superconducting loops with walls thinner than the penetration depth A, finite cur- 
rents are flowing throughout the entire superconductor. The magnetic flux is consequently 
not quantized, but London introduced instead the quantity $' = $ + A/c ^ dr ■ J(r), the 
quantized "fluxoid" . The flux $ is the total flux threading the loop, which already includes 
the current induced flux. A is a phenomenological constant parametrizing the strength of 
the current response of the superconductor to the applied magnetic field; A is related to the 
penetration depth via A = AttX^/c^ through the London equation^. Thin superconducting 
loops therefore react periodically to the continuous variable $. 

It is tempting to relate the $o/2 flux periodicity of superconducting loops to the charge 
2e of the Cooper pairs^ which carry the supercurrent, but the pairing of electrons alone is 
not sufficient to explain the half-integer flux periodicity. A theoretical description of its true 
origin was found independently in 1961 by Byers and Yang^ and by Brenig^ on the basis 
of the BCS theory by realizing that there are two distinct classes of superconducting wave 
functions that are not related by a gauge transformation. An intuitive picture illustrating 
these two types of states is contained in Schrieffer's book on superconductivity—, using the 
energy spectrum of a one-dimensional metallic ring: The first class of superconducting wave 
functions is related to pairing of electrons with angular momenta k and —k and equal energies 
without an applied magnetic field, as schematically shown in Fig. [T] (a). The Cooper pairs 
in this state have a center-of-mass angular momentum q = 0. The wave functions of the 
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FIG. 2: (a) Energy and (b) supercurrent J($) as a function of flux $ for a conventional 

superconducting loop at T = 0. The minima in E{^) correspond to superconducting states with 
different pair momenta q. The screening currents in the superconductor drive the system to the 
closest minimum for each flux value (black points), if the walls of the loop are thicker than A. 

superconducting state for all flux values $, which are integer multiples of $0 ^-nd correspond 
to even pair momenta q, are related to the wave function for $ = by a gauge transformation. 
For a flux value $o/2, pairing occurs between degenerate electrons with angular momenta 
k and —k + 1 [Fig. [T] (b)], and leads to a pair momentum q = 1. The corresponding wave 
function is again related by a gauge transformation to the states for flux values $ which are 
half-integer multiples of $0 and correspond to odd pair momenta. 

The two types of pairing states described above are qualitatively different. For the $o/2 
periodicity, it is further required that the two types of states are degenerate. Byers and 
Yang as well as Brenig showed that this is indeed the case in the thermodynamic limit with 
a continuous density of states. The energy is then determined by a series of intersecting 
parabolae with minima at integer multiples of $0 (corresponding to even pair momenta q) 
and half integer multiples of $0 (corresponding to odd pair momenta q) [Fig. [2] (a)]. If 
the loop is thicker than A, the system locks into the minimum closest to the value of the 
external flux. In flnite systems however, the degeneracy of the even and odd q minima is 
lifted, but their position is flxed by gauge invariance. The flux periodicity in thin loops 
is thus not necessarily $o/2, but the superconducting flux quantum remains $o/2. The 
circulating supercurrent J($) is proportional to dE{^)/d^ and forms a $o/2 periodic saw- 
tooth pattern in the thermodynamic limit as shown in Fig. [2] (b). 
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II. FLUX PERIODICITIES IN CYLINDERS: AN ANALYTIC APPROACH 



For the discussion of the magnetic flux periodicity of rf-wave superconductors we choose 
to bend a discrete two-dimensional N x M square lattice to a cylinder (Fig. [3]) with circum- 
ference Na and height Ma. For two reasons we expect nodal superconductors to support a 
$0 = hc/e rather than a $o/2 periodicity. The first arises from the discrete nature of the 
eigenenergies in a finite system. For the thin cylinder shown in Fig. |3]the mean level spacing 
in the vicinity of the Fermi energy Ep is 6f oc 1/{NM); in s-wave superconductors with an 
order parameter A ^ 6f, 6f matters little. For superconducting states with gap nodes, the 
situation is different. In rf-wave superconductors with an order parameter oc k'^ — k'^, the 
nodal states closest to Ep have to fulfill the condition kz = k^, thus there are fewer possible 
eigenstates and 6f oc 

The second reason is that for gapless superconductors with a finite density of states close 
to E-p, the occupation probabilities of these states change with flux. The flux dependence of 
the occupation enhances the difference of current matrix elements for integer and half-integer 
flux valuesi^""— . This effect is best understood in terms of the spatial extent of a Cooper 
pair. In s-wave superconductors, the occupation probability remains constant for all $, if the 
diameter of the cylinder is larger than the coherence length ^o- If this condition is fulfilled, 




FIG. 3: As a model system we use a thin-wall cylinder constructed from a two-dimensional discrete 
lattice. The interior of the cylinder is threaded by a magnetic flux we assume that the flux does 
not penetrate into the cylinder wall. $ can be chosen arbitrarily, since quantization applies to the 
fluxoid and not the flux itself. 
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the constituents of a Cooper pair cannot circulate separately, i.e. the pair does not feel 
the multiply connected geometry of the cylinder. But for nodal superconducting states, the 
lengthscale which characterizes their coherence, diverges in the nodal directions and there are 
always Cooper pairs which extend around the circumference of the cylinder. Therefore nodal 
superconductors have no characteristic length scale above which the superconducting state 
is unaffected by the geometry of the system. These two combined effects are investigated 
on the basis of an analytic model in Sec. IIIBI 

A. Superconductivity in a flux-threaded cylinder 

The properties of a finite-size multiply connected superconductor depend sensitively on 
the discrete energy spectrum in the normal state. On the N x M square lattice, the flux 
values where levels cross have a high degeneracy for special ratios N/M; for = M, the 
degree of degeneracy is M. For the latter case, the differences between the spectrum for 
integer and half- integer flux values are most pronounced. For N = M ± 1, the spectrum is 
almost •I'o/S-periodic. Away from these special choices of and M, the degeneracies are 
lifted, indicated by the blue shaded patches in Fig. |H The size of the normal persistent 
current circulating around the cylinder is controlled by the change of the density of states 
near Ep upon increasing </> = $/$o- Since normal persistent currents in clean metallic 
rings are typically $o periodio^'^^, we will choose N = M and a half-filled system with the 
chemical potential /i = for our model study, where the $o periodicity of the spectrum 
is most clearly established. Whenever an energy level crosses E-p with increasing flux, the 
current reverses its sign. The current is $o-periodic for even A^ and either paramagnetic 
or diamagnetic in the vicinity of = 0. For odd A^, the current is <l>o/2-periodic. This 
lattice-size dependence persists also in rings with electron-electron interactions^^"— or in 
mesoscopic superconducting islands^. 

We choose in the following A^ and M even, which leads to a normal state spectrum of the 
type shown in Fig. HI This is not an obvious choice, but we will see in chapter 3 that one 
obtains this type of spectrum also for a square loop to which we will compare the results 
obtained for the cylinder geometry. 

The starting point for our analysiss is the BCS theory for a flux threaded cylinder with 
circumference Na = 27TRa and height Ma, where R is the dimensionless radius of the 
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FIG. 4: The energy spectrum of a cylinder in the normal state depends on the numbers and M, 
which parametrize the circumference and the height of the cylinder—. The black lines represent 
the energy levels for a one-dimensional ring with M = 1 and N an integer, where level crossings 
occur for integer values of </> = <&/^o- is the maximum Doppler shift for (j) = 1/2 (see Sec. Ill Bp . 
For Af ^ 1, the levels split and form a quasi continuous spectrum that depends on the ratio N/M 
(blue patches). 

cylinder and a the lattice constant. The pairing Hamiltonian is given by 

H = J] ek(0)cLck. + q)cktc-k+q; + A(k, q)cLk+q44t] ' (2) 

k,s 



where k = (k^.kz) with = n/R and n G {—N/2 + 1, . . . , N/2}. The open bound- 
ary conditions in the z-direction along the axis of the cylinder allow for even-parity so- 
lutions with kz = {2m f. — 1)tt/M and odd-parity solutions with k^ = 27rmo/M, where 
me,mo G {1, . . . , M/2}. The operators cj^.^ and Cks create and annihilate electrons with an- 
gular momentum hk^/a and momentum hkz/a in z direction. For convenience, we choose k^, 
kz G [0,2-7r]. The eigenenergies of free electrons moving on a discrete lattice on the surface 
of the flux threaded cylinder have the form 



ek(0) = -2t cos yk^ ~ ^ 
For i? ^ 1, e]^{(j)) is expanded to linear order in (p/R; 

ek(0) - ek(0) ^ -2t^ sin k^ 



/i. 



(3) 



(4) 



is commonly called the Doppler shift. 

The superconducting order parameter in the pairing Hamiltonian ([2]) is defined through 



A(k, q) = Ag(0)^(k - q/2) = V{\^, k', q) (ck'+q^Ck't - C-k'+qfCk'i)' 



(5) 



where ^(k, k', q) is the pairing interaction. Here we choose a rf-wave interaction in separable 
form: V{'k,k',q) = Vig^k - q/2)^(k' - q/2) with ^(k) = cos(A;^) - cos(A;^); Vi is the 
pairing interaction strength^^. The order parameter A(k, q) represents spin-singlet Cooper 
pairs with pair momentum ^q/a. On the cylinder, the coherent motion of the Cooper 
pairs is possible only in the azimuthal direction, therefore q = (g/i?, 0) with q G {—N/2 + 
1, . . . , N/2}. The quantum number q is obtained from minimizing the free energy. The 
(/)-dependence of Aq(</)) enters through the self-consistency condition and has been discussed 
extensively in^^ anci^ for s-wave pairing, where g(k) = const. As verified numerically, Aq(0) 
varies only little with (p, and we start our analytic calculation with a (p and q independent 
order parameter A(k, q) = A(k) and Ag(0) = A. As in our preceding work—, we take q = 
fioor(2</) + 1/2) in a first step. Since the Hamiltonian is invariant under the simultaneous 
transformation — )■ ± 1 and g — j- g ± 2, it is sufficient to consider g = or 1 and the 
corresponding flux sectors —1/4 < < 1/4 and l/4<0<3/4, respectively. 

The diagonalization of the Hamiltonian ([2]) leads to the quasiparticle dispersion 

(k, q, 0) = '^^^^ - '^-^^^^^^ ± VA2(k)+e2(k,q,0), (6) 

with e(k, q, 0) = [6^(0) -|- e_k+q(0)]/2. Expanding _E±(k, q, 0) to linear order in both (f)/R 
and q/R gives 

E±(k, q, 0) ^ -e,(k) ± y/A2(k) + [ek(0)-/,(k)]^ (7) 

where 

^qi^) = ~ — ^-^2tsink^ and /^(k) = ^ sin /c^. (8) 

In the normal state A = 0, the additive combination of eq(k) and /|j(k) leads to the q- 
independent dispersion ([3]). For A > 0, the dispersion ([7]) differs for even and odd q, 
except for special ratios of and M, as discussed above. This difference is crucial for 
nodal superconductors: The condition ~ kz for levels close to Ep causes a level spacing 
6f ~ 2/i(k^) for small A, where hp is the Fermi momentum. For A^ and M even and q = 0, 
the degenerate energy level at -E = Ep = splits into M levels for increasing A, which 
spread between —A and A. For q = 1, the degenerate levels closest to Ep are located at 
E = ±|/i(ki?)|, thus a gap of 2/i(ki7) remains in the superconducting spectrum. If A^ and 
M are odd, the spectra for even and odd q are interchanged, and if either A^ or M is odd, 
the spectrum is a superposition. 
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The gauge invariant circulating supercurrent is given by 



(9) 



where = de]^{(f)) / d{Rk^p) is the group velocity of the single-particle state with eigenenergy 
ek(0). The spin independent occupation probability of this state is 



n,(k) = (cLcks) = u\k, q, 4>)f{E4k, q, 0)) + v\k, q, 0)/(i?_(k, q, 0)) 



with the Fermi function f{E) and the Bogoliubov amplitudes 
e(k, q,, 



M (k, q, 



E(k,q,, 



+ 1 



and 



V (k,q, 



e(k, q,' 
E(k,q, 



- 1 



(10) 



From Eqs. (Q and (fTOl) . the supercurrent in the cylinder is obtained by evaluating the sum 
either numerically or from the approximative analytic solution in Sec. IIIBt which allows 
insight into the origin of the ^o-periodicity in nodal superconductors. First, the analytic 
solution, which was introduced in Ref.— , is reviewed. 



B. Analytic solution and qualitative discussion 

An analytic evaluation of the supercurrent is possible in the thermodynamic limit where 
the sum over discrete eigenstates is replaced by an integral. For a multiply connected 
geometry, this limit is not properly defined because the supercurrent or the Doppler shift 
vanish in the limit i? — ?■ oo. Care is needed to modify the limiting procedure in a suitable 
way to access the limit of a large but non-infinite radius of the cylinder—. In this limit it is 
mandatory to consider the supercurrent density j(0) = J(0)/M rather than the supercurrent 
J(0). In this scheme, we treat the density of states as a continuous function in any energy 
range where the level spacing is oc 1/NM, but we keep the finite energy gap of width 
2/g(ki?) (X 1/R (X 1/N around Ep in the odd-g sectors. For the tight-binding dispersion 
in Eq. ([3]), the density of states is a complete elliptic integral of the first kind. For the 
purpose of an analytic calculation, a quadratic dispersion with a constant density of states 
is therefore a more suitable starting point. We use the expanded form of Eq. ([3]): 



ek(0) = t 

where fj,' = fj, + it. 



2 

kip ^ ~\~ k.. 



- /i', (12) 
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Some algebraic steps are needed to rearrange the sum in Eq. ([9j) suitably to convert it into 
an integral. For finite 0, ek(0) 7^ e_k(0), and consequently the sum has to be decomposed 
into contributions with > and < 0. We therefore take k^ > and decompose as 

v±k = ^ (^±k^ ~ ^) " '"'^^^'^ ^ ^^^^ 
into a diamagnetic contribution Wd(k) = ~2t(j)/ E? and a paramagnetic contribution fp(k) = 

In a continuous energy integration, the Doppler shift is noticeable only in the vicinity of 
Ey- On the Fermi surface k^p and kz are related by 



k^Ak.) = ^j-kl (14) 

We therefore approximate eg(k) and lq{k) by Cqikz) ~ 2t(0 — q/2)k^pAkz) / R and /^(fcz) ~ 
tqkip^F{kz)/ R, respectively. The eigenenergies ([7]) near i^p are thereby rewritten as 



E+{±k^, kz, q, 0) = Teg(A:,) + V + (ek(0) T /.(A:.)) 



E_(±A;^, kz, q, 0) = Teg(A;,) - ^^k + i^M T k{kz)Y (15) 

The supercurrent 7(0) in Eq. (j9]) is now evaluated by an integral over k^p and kz, which is 
decomposed into an integral over the normal state energy e and an angular variable 9. Within 
this scheme the density of states becomes gapless in the limit M — )■ 00 for g = 0, although 

is kept finite. For q = 1 instead, a /c^-dependent gap 2\li{kz)\ remains. Thus we replace 
ek(0) =F \lq{kz)\ by the continuous quantity e ± \lq{EY,6)\ where we use the parametrization 



kp \ Ik cos 6 




(16) 



kz j \ ksinO I 

with e = tk"^ — fi'. The energy integral extends over the whole tight-binding band width 
with Ep = in the center of the band. Correspondingly, we integrate from —n' to /x'. 
Furthermore, the Doppler shift is parametrized for e ~ Ep as 

eg{9) = ^~'l^'^ 2t^/J7Jicos9 and lq{9) = %^/J7Ji cos 9, (17) 

R R 

where the function lq{9) is positive for |^| < 7r/2. The supercurrent thus becomes 

^^k^ks(q) + ^ fk?^ks(q) 

2^^ [ ^ d9 r de[ng4e,9)v+{e,9)+nq4e,9)v.{e,9)], (18) 
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where nq±{e,9) = 'n'±k(e,6») (q) and v±{e,9) = t'±k(e,6i)- The constant density of states in the 
normal state is A/" = R/Airt. We collect the terms proportional to Vd{^i 0) = —2t(f)/R'^ into a 
diamagnetic current contribution jd and those proportional to Vp{e, 6) = 2tk^^p{e, 9)1 R into 
a paramagnetic contribution jp. Using f{—E) = 1 — f{E), and jp become 

=4Af^r^\er de^,(e,^)^=^=[/(E + e,(e))-/(-E + e,(^))],(19) 
ip(g,0) =AMir'\er dev,{e,e)[f{-E-e,{e))-f{-E + e,m. (20) 

'T'J-7t/2 Jlq{0) 

Here, the integration is over positive e only, and the lower boundary of the energy integration 
is controlled by lq{9). In Eq. ( fT9l) we used the abbreviations A = A(^) and E = E{e,9) = 
A'^(6) + e^. The current jd is diamagnetic in the even-g flux sectors and paramagnetic 
in the odd-g sectors. For even q, it is equivalent to the diamagnetic current obtained from 
the London equations^^i^. The current jp has always the reverse sign of ja and is related 
to the quasiparticle current as shown below. To analyze the flux dependent properties of 
the spectra and the current in the even-g and odd-g sectors, we explicitly distinguish s-wave 
pairing and d-wave pairing with nodes in the gap function. 



1. s-wave pairing symmetry 

For s-wave pairing, A(e,^) = A is constant. Therefore, if we assume that A > eq{6) for 
all 6, the lower energy integration boundary in Eqs. f|T9l) and fl20l) is A. Thus j(0) = Jd+Jp is 
equal in both the even-g and the odd-g flux sectors and the flux periodicity is $0/2. However, 
if A < max0eg(6'), different calculational steps have to be followed in the evaluation of 
Eq. ([9]), the results of which have been presented 



With e = \/E'^ — A^, Eqs. ( IT9l) and (l20l) transform into integrals over E with de = 
D,{E) dE, where 

de ( E (E^ - A^)'^/^ for E>A 
DAE) = ^={ ^ - (21) 

^ ^ 5^ \o for E<A ^ ^ 

is the density of states for s-wave pairing. This leads to 

jd = ^\ f^'^ de £ dEv, [y/W^, e) [f{E + eq{e)) - f{-E + 6^(6))] , (22) 
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f7r/2 

At T = 0, we obtain 



-W^ / d9cos9 dEDs{E) + 0(-) . (25) 

- V t J-tt/2 J a 



-7r/2 

StATe^ /77 /-"/^ ^^^^^^ 

The current jd becomes independent of the superconducting density of states. Its size is 
proportional to i?F, as long as /i' ^ A holds. 

If A > eq{6) for all values of 6, then jp = and the supercurrent = jd is diamagnetic. 
For T > 0, jd decreases slightly. The current jp increases with increasing T and reaches its 
maximum value at Tc. For finite temperatures jp is referred to as the quasiparticle current. 
The supercurrent is always the sum of the diamagnetic current jd and the quasiparticle 
current jp, and therefore decreases with increasing temperature and vanishes at Tc— . The 
quasiparticle current has the same flux periodicity as the supercurrent, even though it is 
carried by quasiparticle excitations. In the normal state (A = 0), 

(26) 

which cancels jd exactly in the limit^ M — > 00. 



2. Unconventional pairing with gap nodes 

Equation (124|) for jd is valid also for unconventional order parameter symmetries. Phys- 
ically, jd reflects the difference in the density of states of quasiparticle states with orbital 
magnetic moments parallel and anti-parallel to the external magnetic field. The former 
states are Doppler shifted to lower energies, whereas the latter are Doppler shifted to higher 
energies. This is schematically shown in Fig. [5] for d-wave pairing (c.f.— ). In this picture, 
jd is proportional to the difference between the area beneath the red and and blue curves 
representing the density of states arising from i?_(ih|k|, q, </>) < 0. Therefore we approximate 
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a) 



Uj, 




b) 



Uj. 




FIG. 5: Scheme for the density of states of a d-wave superconductor for (p = 1/4, where Cq = Zi/22i. 
The center-of-mass angular momentum hq/aR of the Cooper pairs is (a) q = and (b) q = 1. The 
energies are Doppler shifted to higher (red) or lower energies (blue). This results in a double-peak 
structure; for q = the upper and lower band overlap in the region — eo < E < ecr^ and states in 
the upper band become partially occupied. For q = 1 there is a gap h of the size of the maximum 
Doppler shift at i;^> = 1 /4. The black line represents the density of states (a) for </> = and (b) for 
= 1/2. 



jd for A(^) < /i' = + 4t by 



3d 



-2/i'- 



q/2 



(27) 



h R ' 

as given in equation with A = 0. On the other hand, jp is represented by the occupied 
quasiparticle states in the overlap region of -E+(k, q, </>) and -E-(k, q, (p) with width 2eq(kF)- 
It therefore strongly depends on the density of states in the vicinity of Ep. In Fig. |5] (a), 
which refers to even g, the current jp is determined by the small triangular patch where the 
upper and lower bands overlap. For odd q, the two bands do not overlap, therefore jp = 0. 

We will now analyze such a scenario for d-wave pairing with an order parameter Ak = 
A(fc^ — k"^) ^ A cos 26'. Again, we assume A > eq{6) for all 6; then the integral in Eq. (120|) 
contains only the nodal states closest to Ep, for which the (i-wave symmetry demands k^p ^ kz- 
Jointly with Eq. iHM this condition fixes the Doppler shift at Ep to the k-independent value 



- q/2)y/2tJ7/R and Iq = {q/ R)^/iJ7j2. With the density of states 

Dd{E) 



VE^-A^cos^ 29 
Eq. (l20l) for the paramagnetic current jp at T = takes the form 



(28) 



e f" 



Tv/2 



■k/2 



d9D,{E) 



R 



2t e + fi' . 



sin 9. 



(29) 
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FIG. 6: The supercurrent density = jd + ip in a thin d-wave cyhnder as a function of flux 
(arbitrary units}^. Shown is the result of the analytic model calculation [Eq. ()33p ] for the 
characteristic value h = 0.4. For —1/4 < (j) < 1/4, where q = 0, the current is reduced by a 
contribution proportional to cjP' , whereas it is linear in (j) otherwise. This leads to an overall flux 
periodicity of <I>o- 

In the odd-g flux sectors, Iq > Cq for all values of 0, therefore jp = 0. In the g = sector, 
/g = and 



Tv/2 



tt/2 



- A2 cos2 26 7lh\ t Jo A 
TThAy t TtA^^ h V i? 



(30) 



where the same approximations as in the s-wave case are applied. The dominant contribution 
to the integral over 6 originates from the nodal parts (see e.g.—). 
In the even-g sectors, the total current j(0) = jd + jp becomes 



vrA R 

which results in the ratio of the two current components 



(31) 



Jp 
3d 



ttA R 



hcj). 



(32) 



In the odd-g flux sectors jp = and the supercurrent is j(0) = jd- j(0) is consequently 
$0 periodic; within one flux period from —1/2 to 1/2 we represent it as 



'Rh 



+ 1/2 for -1/2 < (p< -1/4, 
;i - h(j)) for -1/4 < (j)< 1/4, 
- 1/2 for 1/4 < < 1/2, 



(33) 
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(c.f. Fig. [6]). The difference of the supercurrent in the even-g and odd-g ffux sectors is 
represented best by the Fourier components j„ = J^y2^^ ^i^)^'^''^"^'^- -^^^ ^^^^ (^i) 
the second Fourier component (j^) we obtain 

/i' e ini — b 



3i 



IT 



and 



32 = -2'- 



(34) 



Rh 167r3 Rh 167r2 ' 

To leading order in 1/R, the ratio of the $o a-^d the $o/2 Fourier component is therefore 

2t 



hm 

/i-)-0 



0.07 



(35) 



47r2 Ai? ' 

and scales with the inverse ring diameter. This l/i?-law is the direct consequence of the 
(i-wave density of states Dd{E) oc E. Using Eq. ( l35i) to estimate this ratio for a mesoscopic 
cylinder with a circumference Ra = 2600a ~ 1 /xm and a ratio A/t = 0.01, we obtain 
Ji/J2 ~ 0.03. 



C. Further aspects 

We have shown that in rings of unconventional superconductors with gap nodes, there is 
a paramagnetic, quasiparticle-like contribution jp > to the supercurrent at T = 0. This 
current is generated by the fiux-induced reoccupation of nodal quasiparticle states slightly 
below and above Ep. Formally a coherence length /ifir/A(k, q) > 2R can be ascribed to 
these reoccupied states, which are therefore affected by the symmetry of the system. If 
the normal state energy spectrum has a flux periodicity of $0; than the superconducting 
spectrum is $o periodic, too. The normal state spectrum of a cylinder with a discrete lattice 
strongly depends on the number of lattice sites. This problem is characteristic for rotationally 
symmetric systems and is much less pronounced in geometries with lower symmetry, such 
as the square frame discussed in Sec. 3. In the latter system impurities do not change the 
spectrum qualitatively. For modelling an experimental arrangement a square loop geometry 
is therefore preferable. 

The $0 periodicity is best visible in the current component jp at T = 0. For d-wave- 
pairing jp oc l/R"^, and the $o periodic Fourier component decays like the inverse radius of 
the cylinder, relative to the $o/2 periodic Fourier component. The lack of a characteristic 
length scale in nodal superconductors, such as the coherence length for s-wave pairing, 
generates this algebraic decay. Although jp is larger for small A, it almost vanishes close to 
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Tc, if A ^ Sp, and variations of Tc with flux, as in the Little-Parks experimental^, do not 
differ for s- and d-wave superconductors. 

III. FLUX PERIODICITY IN SQUARE FRAMES: BOGOLIUBOV DE GENNES 
APPROACH 

So far we have presented the principles of the crossover from $o to $o/2 flux periodicity 
in conventional and unconventional superconductors and the mechanisms that leads to the 
persistence of $o periodicity in large loops of nodal superconductors. Now we present an 
alternative approach in real space via the Bogoliubov - de Gennes equations, which we 
introduce in Sec. Illl A[ The information we obtain from this technique is complementary 
to Sec. ini where we followed the momentum-space formulation. The latter proved useful 
to understand the physical concepts and to describe large systems. The price paid was 
the restriction to highly symmetric systems with intriguing energy spectra in the normal 
state. This raises the question whether the $o periodicity is detectable in realistic setups, 
or whether it is rather an artifact of the high degeneracy of energy levels in clean and highly 
symmetric systems? On the other hand, the Bogoliubov - de Gennes equations in real space 
allow to determine the spectrum of "natural" system geometries with reduced symmetry or 
systems containing lattice defects, impurities, magnetic fields or correlations in real space. 
Limitations of computational power, however, restrict the system size, and therefore the 
particular effects introduced by discreteness are unavoidably present. 

The combination of momentum- and real-space methods can provide answers to the 
questions above. In the following, we first discuss the multi-channel loop for a square 
lattice: a square frame, as shown in Fig. U\ with a square hole at the center, threaded by 
a magnetic flux $. We use this system in Sec. IIII Bl to study the flux periodicity in clean 
symmetric square frames; a part of this section is contained in^. In Sec. IIVI we investigate 
different Josephson junction devices that respond periodically to magnetic fields. Junctions 
are modeled in real space by inserting potential barriers. In this context, we investigate also 
the effect of impurities and lattice defects on the energy spectrum of the square frame. 
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A. The Bogoliubov — de Gennes equations 



The Hamiltonian which we use in the following section has the form 



(36) 

where c|j,, are creation and annihilation operators for an electron on lattice site i with 
spin s, and /i is the chemical potential. The sum runs over all lattice sites and the sum 
is restricted to nearest-neighbor sites i and j only, and tij = te'^^^ with the hopping 
amplitude t and the Peierls phase factor 



^., = ^^'dr.A(r). (37) 

Additionally, we include an impurity term consisting of potential scatterers with repulsive 
potentials Ui > 0, which we align to model tunnel junctions. A Hamiltonian of the form ( 136|) 
has often been used before for the numeric investigation of vortices in d-wave superconductors 
and the technique is described in detail in a number of articles^"—. 

In the Hamiltonian Eq. fl36|) two types of spin-singlet pairing are included. The on- 
site order parameter Aj represents conventional s-wave pairing originating from an on-site 
interaction. The order parameter originates from a nearest-neighbor interaction between 
the sites i and j. They are defined through 

Ai = Vo(c4Ci|) and Aij = y [(cj^Qt) - (Q|Cjt)] • (38) 

with the interaction strengths Vq and Vi. To diagonalize the Hamiltonian fl5^ we use the 
Bogoliubov transformation 



At 



n n 

where the coefficients Un, and Vn, are obtained from the eigenvalue equation 




E^l^]. (40) 



'n 



The operators t and A act on the vectors u,„ and v,„ as 

= ^UjUnj + {Ui - n)un, and Avn, = Ait;„, + ^ AijVnj, (41) 
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where j labels the nearest-neighbor sites of site i. Inserting the transformation (139|) into 
Eq. fl38|) leads to the self-consistency conditions 



Ai = VoY^ Mn,< tanh (^j , Aij = y J] [^^^nj + Wnj<] tanh (^j . (42) 

Equations ( H2|) together with Eq. ( HOj) represent the Bogoliubov - de Gennes equations. 

The bond order parameters Aj^ can be projected onto a d-wave component and an ex- 
tended s-wave component defined as 

At = ^ [A,,i+£e'^-+* + Ai,i_£e'^--^ + A,,,+^e''^-+^ + Ai,,_^e*'^-'-^] , (43) 

Af = 1 [A,,,+£e^'^-+* + A,,,_£e^'^'-* - A,,,+ge''^-+^ - A,,,_ge''^--] . (44) 

In a uniform system with nearest-neighbor pairing interaction only, the self-consistency 
Eq. ( H2l) selects a pure d-wave superconducting state, i.e. A| = 0. Impurities, potentials or 
boundaries generate an extended s-wave contribution A| > 0^. The expectation value of 
the current Jij (cf.— ) from site i to j is given by 

J,, = -8t $0 V Im {unjule-"^^^) f{E^). (45) 



B. Flux periodicity in square frames 

The Bogoliubov - de Gennes equations introduced above are now applied to the square 
frame geometry shown in Fig. [71 consisting of a discrete N x N lattice with a centered Lx L 
square hole threaded by a magnetic flux (f), where </> = $/$o- The external magnetic field B 
threading the hole is supposed not to penetrate into the frame, and we restrict it to the center 
of the hole. B is generated by a vector potential of the form A(r) = 27r0/|rp(?/, —x, 0). 

In the normal state the Bogoliubov - de Gennes equations reduce to the discrete Laplace 
equation. While the low-energy states do not differ much from free plane waves, the higher- 
energy states near Ep on the square frame develop some peculiar, frame-specific features. 
The wavelength of a state near Ep is close to two lattice constants, therefore the probability 
density divides into two sublattices. In the square frame, structures on different sublattices 
can overlap, which results in the characteristic real-space density profiles which persist in 
the nodal states of a d-wave superconductor. Figure M shows two such examples. 
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FIG. 7: Illustration of a square loop threaded by a magnetic flux. For the investigation of the flux 
periodicity of d-wave superconductors using the Bogoliubov - de Gennes equations in real space, 
we use a discrete square lattice with open boundary conditions and a square hole in the center the 
frame, which is pierced by the magnetic flux 

The characterization of the superconducting solutions of the Bogoliubov - de Gennes 
equations in the square frame is analogous to those on the cylinder in the momentum space 
analysis. The absolute value of the d-wave order parameter \Af\ is shown in Fig. [9] (a) for 
= 0. The open boundary conditions cause a decrease on the boundaries and are respon- 
sible for Friedel oscillations visible along the diagonal. In multiply connected geometries, 
the Bogoliubov - de Gennes equations generally allow for solutions where acquires a 
phase gradient such that the phase difference on a closed path around the hole is 27rg with 
integer q. As in Sees. [T] and [Tll this phase winding number q represents the center-of-mass 
motion of a Cooper pair, although it cannot be identified with the angular momentum in the 
square geometry. The different numerical solutions are obtained by choosing appropriate 
initial values for the phase of A^, and the phases of the self-consistent results are shown in 
Figs. M (b), (c) and (d) for g = 1, 2 and 3 and flux values = 1/2, 1 and 3/2, respectively. 

To assess the E{(j)) and the current J{4>), the evolution of the eigenenergies with magnetic 
flux has to be calculated first. The eigenstates with energies below Ep form the ground-state 
condensate (Fig. ITOl) . Here we discuss only flux values </> between and 1/2, because all 
quantities are either symmetric or antisymmetric with respect to flux reversal — )■ —(f). The 
spectrum for a square frame with = 40 and L = 14 is shown in Fig. [10] for half flUing, i.e., 
/i = 0. Because the number of lattice sites on straight paths around the hole is a multiple 
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of four in a square frame, the spectrum is almost identical to the one for a cylinder with an 
even number of lattice sites and with the same number iV — L = 26 of transverse channels 
(compare to Fig. 6 in Ref.— ). For the square frame, the energy levels do not actually cross 
Ep, because the lack of rotational symmetry leads to hybridization of the levels and level 
repulsion. Nevertheless, the same clearly distinct flux regimes are found: the flux intervals 
between and 1/4 and from 1/4 to 1/2 (in units of $o)- 

Up to ~ 1/4 the current J{4>) generates a magnetic field which tends to reduce the 
applied field by a continuous shift of the eigenenergies in the condensate. At = 0, pairs of 
states with opposite circulation compensate their respective currents, thus J = 0. The well 
separated states at = in Fig. [TU] are the states in the vicinity of the nodes of the d-wave 
superconductor. Away from Ep, the density of states increases towards the states near the 
maximum energy gap A that provide most of the condensation energy. For > 0, the 
energy of the states with orbital magnetic moment anti-parallel (parallel) to the magnetic 
field is increased (decreased). Correspondingly the supercurrent, which is carried by these 
states, depends on the details of level crossings and avoidings. The main contribution to the 
supercurrent arises from the occupied levels closest to Ep, because the contributions from 




FIG. 8: Real-space representations of a square loop with a typical electronic probability density 
l^p. We show two eigenstates of the d-wave pairing Hamiltonian with slightly different energies 
in the gap region, calculated for a square- loop with 80x80 lattice sites and a pairing interaction 
Vi = 0.3t. The hole in the center has a size of 28x28 unit cells. To enhance the contrast of the 
complicated pattern, the special color code shown on the right is used and the discrete lattice 
points are smoothly interpolated. 



20 



40 



40 




FIG. 9: (a) Absolute value of the d-wave order parameter (in units of t) in a 40 x 40 square 
frame with a 14 x 14 hole at the center for q = 0, (p = 0, and Vi = 0.3t. For this interaction 
strength, the suppression of Af for (j) ^ is small and not visible in this plot. The phase of Af is 
shown for winding numbers g = 1,2,3 in (b), (c) and (d), respectively. 

the lower- lying states tend to cancel in adjacent pairs. 

As the highest occupied state shifts with increasing flux to lower energies, the current in 
the square loop first increases for small (Fig. fTTl) . then decreases when the highest occupied 
level with an orbital moment opposite to the applied magnetic field starts to dominate. With 
increasing flux this state approaches E-p. A current-carrying state in the vicinity of the nodes 
is replaced upon a slight increase of by a state of opposite current direction. The states of 
the condensate are thereby continuously changing near the extrapolated crossing points. As 
a consequence, the energy "parabola" centered at zero flux is different from the ground-state 
energy parabola centered at </> = 1/2 [Fig. [TT] (a)]. The deviation from a parabolic shape 
near zero flux is due to the evolution of the near-nodal states; the vertical offset of the energy 
minima at (p = n results mostly from the flux dependence of the states near the maximum 
value of the anisotropic gap. 
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FIG. 10: Energy spectrum for a d-wave superconductor on a square frame. The eigenenergies in 
the gap regfon are shown for a square 40x40 foop with a 14x14 hole and pair interaction Vi = 0.3t 
as a function of flux (j) (in units of <^o)- The energies are given in units of the superconducting 
order parameter A at = (A ~ 0.22t). The superconducting condensate consists of the states 
below Ep = 0. Reconstruction of the condensate takes place near (p = ib(n + l)/2, where the 
eigenenergies jump abruptly (after Ref.— . 

For flux values near </> = 1/4 the condensate reconstructs. The superconducting state 
beyond 1/4 belongs to the class of wave functions introduced by Byers and Yangi^ in which, 
for a circular geometry, each pair acquires a center-of-mass angular momentum hr^. Re- 
markably, in the flux interval from near 1/4 to 1/2, a full energy gap exists also for d-wave 
superconductors (Fig. [TOl) . Here the circulating current enhances the magnetic fleld; the 
paramagnetic orbital moment of the current is parallel to the fleld. The resulting energy 
gain is responsible for the fleld-induced energy gap. This reconstruction of the condensate 
is the origin of the $o periodicity in energy and current. 

These calculations show that a rf-wave superconducting loop in a square geometry has 
almost identical properties to a flux threaded cylinder. This is remarkable, because on a 
closed path in the square frame, the phase of the ci-wave order parameter Af rotates by 
27r, whereas in the cylinder, the order parameter rotates with the lattice. Therefore, while 
changes in the geometry and the number of transverse channels modify the spectrum and 
the J{4>) characteristics in detail, they do not eliminate the $o periodic component. The 
reduction of the symmetry, here to the four-fold rotational symmetry of the square frame, 
stabilizes the spectrum compared to the cylinder geometry. 
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FIG. 11: Flux dependence of energy and current for the square frame. Energy [E{(j)) — E{0)]/E{0) 
(a) and circulating current J{(p) (b) for a square 40x40 loop with a 14x14 hole and pair interaction 
Vi = 0.3t. J(0) is given in units of t/^o = 6 x 10~^A for the choice of t = 250 meV. The 
condensate states with even and odd winding number q are clearly distinct, which is reflected, e.g., 
in the deformation of the q = 0-parabola. The overall (j) periodicity for E{(f)) and J{(f>) is $0 (from 
Ref^). 

IV. FLUX PERIODICITY OF JOSEPHSON JUNCTIONS 

All energy levels are Doppler shifted in current carrying systems, not only in flux threaded 
loops but also in wires or at the surface of bulk superconductors. In the latter systems, the 
phase gradient of the superconducting order parameter typically does not reach the value nec- 
essary to drive the superconductor into a finite-momentum pairing state with q ^ 0, which is 
why the influence of finite momentum pairing on the flux periodicity has not been discussed 
in the literature until recently. An exception are systems with strong inhomogeneities of the 
order parameter, which act as Josephon junctions. The phase gradient accumulates at the 
junctions and they behave periodically with the phase gradient, as described by the Joseph- 
son relation. From what has been discussed for the flux periodicity in multiply connected 
geometries, it appears natural that the Doppler shift of nodal states might also influence 
the periodicity of Josephson junctions. 

A Josephson junction is intrinsically a more complicated system than a superconducting 
loop. Several parameters are needed to characterize the junction as well as the superconduct- 
ing states on each side of the junction. Most junctions can be classified either as transparent 
or as tunnel junctions, regardless whether they consist of a geometrical constriction, a po- 
tential barrier, or a normal metal bridge. This classification is closely related to the Doppler 
shift of single energy levels in the system, as will be explained below. In the following we will 
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therefore discuss the Josephson relations in both the tunnehng and the transparent regimes. 
A. Current-phase relation 

The current-phase relation, which expresses the supercurrent J over a Josephson junction 
as a function of the phase difference 5ip of the order parameters on both sides of the junction 
is: 

J = JcSin(5(/9). (46) 

Jc is the critical current over the junction, above which the zero voltage state breaks down. 
This relation was predicted by Josephson in 1963^ and can be directly derived from a 
Ginzburg-Landau description^^. For transparent junctions, in Eq. fH6l) distorts into a 

saw-tooth pattern similar to the current-flux relation in superconducting loops^. It is crucial 
to realize that the phase gradient of the order parameter is twice that of the superconducting 
wave function. If the phase difference of the order parameter on both sides of the junction 
is 5if, then the phase difference of the wave function is 5(^/2. Because the wave function 
of the system has to be 27r-periodic, the periodicity of the energy spectrum and the order 
parameter of a finite system is 4-7r. The current contributions from all energy levels add up 
to a 27r periodic supercurrent only in the thermodynamic limit. In this section we analyze 
whether the Doppler shift of the energy levels leads to the same doubling of the periodicity 
in 5ip of a junction as it does for the flux periodicity of loops. While for the tunneling regime 
we rely on a simple linear-junction model, we will analyze transparent junctions by inserting 
a Josephson junction into a square frame. This has the advantage of a remarkable stability 
of the energy spectrum against the insertion of impurities and lattice defects, as will be seen 
in Sec. IIVA2I 

1. Tunnel junctions 

A simple model of a tunnel junction is a square lattice with sites in x-direction and 
M sites and periodic boundary conditions in ^/-direction. The junction is modeled in the 
tunneling regime by one or two lines of potential scatterers with a repulsive potential U > At 
(Fig. [T2l) . In the absence of a magnetic field, this system is homogeneous in y-direction, and 
the Fourier transformation with respect to the y-coordinate will allow the diagonalization 
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FIG. 12: To model a Josephson junction we choose a discrete square lattice with N sites in 
rr-direction, and M sites in y-direction. The junction itself is modeled by one or two lines of 
potential scatterers (black points) each with a repulsive potential U > 0. 

of larger systems^. 

The Bogoliubov - de Gennes equations are slightly modified in this case: the eigenvalue 
equation fl40l) for nearest-neighbor interaction is now defined through the relations 

iUn , = / UiUn , + (Cfc + Ui )Un i. , Aw„. , = > Aj.,-t;„, , + Afc , , (47) 

Jx Jx 

where eky = —2t cos{ky)—fi and A^^^ = Ay cos{ky). The indices of the eigenvectors u„ and v„ 
are the x-coordinate of the site and the wave number ky in y-direction. The corresponding 
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FIG. 13: Current-phase relation calculated for Josephson junctions in the tunneling regime. Left 
panel: = 18, M = 12, Vi = 0.3t, and U = 4.5f. The sm.{5ip) relation is considerably deformed, 
which is typical for narrow junctions with very few channels. Right panel: N = 20, M = 200, 
Vi = t, and U = 5t. This junction has a sufficiently many channels to exhibit the known current- 
phase relation. The overall sign depends on the choice of the sign in the phase factor of the hopping 
matrix elements. 
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self-consistency equations are 



\i = -y \un- , < + . < 1 tanh ( ^^^) , (48) 

if j = z ± X, and if the bonds are along the y direction 

A. = 1^1 E «n...,<^,^ cos(fc,) tanh . (49) 

The self-consistency equation for the s-wave order parameter Aj with on-site interaction is 
analogous to poj) . but without the factor cos(A;y). 

To induce a finite phase gradient of the order parameter and a supercurrent, we introduce 
a "phase jump" 5ip in the matrix elements tij for hopping from = N — 1 back to = 0, 
and a jump —6lp for the corresponding hopping in the opposite direction. An alternative, 
but physically equivalent choice for the phase of tij is a constant phase factor 6**^'^ with 
ifij = 6ip/N for all hopping processes along the x-direction, which is mathematically identical 
to a cylinder threaded by a flux $ = (hc/e) (p with 27r0 = 6(f/2. In the fully transparent case 
with U = 0, this leads to a homogeneous phase gradient of Af = (Aj^j+£ -|- Aj^j_f)/2 + A^ 
(or Aj, respectively), whereas far in the tunneling regime for U > 4t, the phase of the order 
parameter drops only across the junction. The current across the junction is calculated as in 
Eq. ( H5l) . The results for two typical situations are shown in Fig. [131 The left panel displays 
the current-phase relation of a narrow Josephson junction with a width of M = 12 sites. The 
usual current-phase relation is considerably deformed in this CcLSG, clS IS typical for junctions 
with very few channels^. The exact form of the current-phase relation is characteristic for 
each junction; it depends on the structure of the energy spectrum, which changes strongly 
upon increasing or decreasing the system size or adding impurities. For increasing M, the 
current-phase relation approaches (H6|) . as the level spacing becomes negligible. This is the 
regime of wide junctions, shown in Fig. [13] (right panel), which is well described by the 
Ginzburg-Landau approach. 

Our numerical analysis shows that the Josephson relation fH6l) describes wide junctions 
in the tunneling regime very well; a doubling of the period is not observed, even for d- 
wave superconductors with small antinodal energy gap. The reason for this is twofold: (1) 
Along with the suppression of the critical current Jc across the junction, the Doppler shift 
decreases strongly with increasing repulsive potential U. In the tunneling regime {U > At), 
Jc decreases by a factor > 10'^. Consequently no energy levels (or negligibly few in very 
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FIG. 14: The top 14 energy levels below Ep calculated for a linear Josephson juntion with N = 18, 
M = 12 and Vi = 0.6t as a function of the repulsive potential U on the junction. 

large systems) approach Ep as a function of 6ip and the effects related to a reversal of single 
particle currents are absent. (2) For tunnel junctions, the thermodynamic Ginzburg-Landau 
limit is reached also for (i-wave superconductors, if the density of states close to Ep becomes 
quasi-continuous, in contrast to the flux threaded loop. The deformation of the current- 
phase relation in narrow tunnel junctions is generically not due to levels reaching Ep. The 
deformation is induced, if the total current is carried by very few states, each with a period 
of Aw. The 2 TT- asymmetric terms do not cancel and the critical current is 47r-periodic. 



2. Transparent junctions 

Transparent junctions are more involved than tunnel junctions. One reason for their 
complexity is the strong coupling of the superconducting states on both sides of the junc- 
tion, which does not allow to choose the phases of the corresponding order parameters inde- 
pendently. Consequently the phase difference 6ip is not an adequate variable for describing 
the current across the junction. Another reason is that the energy spectrum in a linear 
junction of the type shown in Fig. [12] changes strongly upon changing microscopic details of 
the system, such as the strength of the repulsive potential U on the impurity sites in this 
case. This is illustrated vividly by Fig. [14] showing the evolution of the highest occupied 
energy levels with increasing U. 

These problems can be resolved by using a square-frame geometry as in Sec. IIIIBI Here 
the Josephson junction is modeled by adding potential scatterers on a line as shown in 
Fig. [ini and the current is driven by a magnetic flux (p threading the frame. For a tunnel 



27 



junction, this would induce a phase jump of 47r0 in the order parameter aross the junction 
and thus a sin(47r0) current-flux relation. In transparent junctions, the jump is smaller and 
vanishes in a clean frame. For this topology the magnetic flux $ = ■ hc/e is related to the 
phase variation of the order parameter 5if by 27r0 = 5ip/2. 

For sufficiently large [/, say U = lOOt, these impurities act as a geometrical constriction. 
Figure [16] shows explicitly that the spectrum of a square frame remains qualitatively invariant 
upon inserting a small number of impurities, even sufficiently strong to block the current 
over the impurity site completely. Figures [16] (b) and (c) show the spectra versus for two 
and four impurities for a 20 x 20 square frame with a 8 x 8 square hole. In the presence of 
impurities, bound states arise at Ep in a (i-wave superconductor— i^, which are nearly flux 
independent. These bound states are easily identified in Figs. [Tn](a) and (b) near the Fermi 
energy. Otherwise, the spectrum in Fig. [16] (b) is very similar to that of the clean frame 
discussed in Sec. IIII Bl (Fig. [TOj) . Clearly visible is the discontinuity of the spectrum where 
the condensate reconstructs to a superconducting state with different winding number q. 
The relevance of g is a characteristic property of transparency and directly connected to a 
discontinuity of the supercurrent (see Fig. [16] (a) for one and two impurities). 

For three to five impurities, the supercurrent is continuous, as is the spectrum shown 
in Fig. [16] (c) for four impurities. Nevertheless, the typical features of the square-frame 
spectrum are still present, in particular the gap in the odd flux regimes and one energy 
level approaching Ep in the even q regime. This level causes the wiggle in the supercurrent 
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FIG. 15: For the description of transparent junctions, we choose the square-frame geometry and 
model the junction with potential scatterers arranged on a Une crossing one side on the frame 
(black points). The current is driven by a magnetic flux threading the frame. 
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FIG. 16: Supercurrent J{(j)) and energy spectrum E{(f)) of a 20 x 20 square frame with an 8 x 8 hole 
containing a Josephson junction. The width of the arms for the frame is six sites. The impurity 
potential is U = WOt. (a) J{(j)) for one (blue), two (turquoise), three (green), four (orange), and 
five (purple) impurity sites. Energy spectrum for (b) two and (c) four impurity sites. 

around (j) = 0; its slope and that of a few others remain almost as steep as in a clean frame, 
which indicates the existence of channels with free current flow. The Doppler shift of nodal 
states is therefore not negligible in the calculation of the supercurrent across transparent 
Josephson junctions, and it may cause appreciable deviations from the sin(47r0) current-flux 
relation even in the case of wide junctions. The sin(47r0) is expected for the thermodynamic 
Ginzburg-Landau limit. 

Finally we note that for five impurities, only one channel through the junction remains, 
which is almost blocked by the bound state. Thus the spectrum becomes nearly flux inde- 
pendent, leading to a junction in the tunneling regime. However, the supercurrent does not 
follow the expected sin(47r0) but rather a sin(27r</)) current-phase relation. This is due to 
the point-contact like character of the junction and the extreme limit of the deformation 
of the current-flux relation as shown in Fig. [13] - similar to the left panel, however with 
J{6ip) ~ -sin(5v9/2). 
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B. Field-threaded junctions 



A magnetic field threading a Josephson junction modifies the phase difference of the order 
parameters of the superconductors on both sides and thus alters the supercurrent. This 
behavior is well understood on the basis of the Ginzburg-Landau approach. The current- 
flux relation of a linear junction that is homogeneous in y-direction has the shape of a 
Fraunhofer diffraction pattern^^, although it deviates from the Fraunhofer form for all other 
junction geometries. Despite these deviations it preserves the characteristic flux periodicity 
of $0/2 for conventional Josephson junctions. The magnetic field dependent critical current 
of Josephson junctions is therefore another key property where the Doppler shift might cause 
a doubling of the flux period. 

Here we use again the linear junction model of Sec. II V A H and fix the phase difference to 
5ip = 7r/2, for which the absolute value of the current across the junction in the tunneling 
regime is largest. In order to introduce a magnetic field threading the junction, we construct 
the junction from single plaquettes with potential scatterers on each of its sites. All plaquette 
/ which belong to the junction are threaded by a magnetic flux generating Peierls phase 
factors (fly We restrict our discussion to a homogeneous field distribution inside the junction, 
(f)^ = (j) for all /, and the repulsive potential on the respective sites is U. In the presence of a 
magnetic field, the system is not homogeneous in y-direction, and we have to diagonalize it 
in real space. This restricts the maximum system size for our analysis. 

1. Current-flux relation of tunnel junctions 

The simplest model of a field-threaded Josephson junction consists of two lines of impurity 
sites as used in Sec. IIV A II (Fig. [T5]). The current-flux relation of such a junction as obtained 
from the Bogoliubov - de Gennes equations is shown in Fig. [TTjfor s- and (i-wave junctions 
with a length of 14 sites and thus 13 plaquettes. Upon first glance, the current-flux relation of 
the s-wave junction [Fig. [17] (a)] appears to be similar to the Fraunhofer pattern known from 
the Ginzburg-Landau approach for linear Josephson junctions^S, as does the current-flux 
relation for the (i-wave junction [Fig. [T7| (b)]. The characteristics are a central peak around 
= with width $0 and side peaks of decreasing height with width $o/2- They display 
the expected global periodicity of 13$0) enforced by gauge invariance, if each plaquette is 
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FIG. 17: Absolute value of the maximum current flowing across a tunnel junction versus the total 
apphed magnetic flux (p (in units of hc/e) obtained from the Bogohubov - de Gennes equations 
solved on a lattice with = 16, M = 14, U = 5t, and 5ip = 7r/2 in units of J(0). (a) s-wave 
pairing with Vq = t. (b) d-wave pairing with Vi = 0.7t. Here (j) is the flux within the junction. 

threaded by an integer multiple of $o- On closer inspection of Fig. [T7] (a) one finds that 
the s-wave junction has one maximum surplus in one period of 13$0) whereas the d-wave 
junction has not. The width of the peaks in Fig. [T7| (a) is therefore slightly smaller than 
the expected value $o- In the following, we explain this effect jointly with an investigation 
of the current-flux relation of inhomogeneous junctions by analyzing the Ginzburg-Landau 
approach for a lattice model. 

We consider a two-dimensional superconductor which is divided by a thin, quasi one-di- 
mensional Josephson junction of width d oriented along the y-direction with <^ A, such 
that screening currents are negligible; A is the London penetration depth. If the junction is 
threaded by a constant magnetic field Bz{x,y) = B^, the supercurrent across the junction 
derived from the Ginzburg-Landau equations is 



J 



dyjc{y) sm{ky), 



(50) 



where k = nBzd/^Q. The critical current density jc{y) is controlled by the microscopic 
structure of the junction. If jc{y) is constant, one obtains the well known Fraunhofer pattern 

(51) 

for the current across the junction; $ is the total magnetic flux through the area of the 
junction. 

On a discrete square lattice with M lattice sites in ^/-direction and an order parameter 







sin(7r<l>/<l>o) 


J(0) 




7r<l>/<l>o 



31 



a) 



\ 












w 













b) 



10 




FIG. 18: Current-flux relation of a Josephson junction as in Fig.[T7]but with M = 11 (10 plaquettes) 
as obtained from the BogoUubov - de Gennes equations (blue) and from the discrete Ginzburg- 
Landau approach (turquoise), (a) s-wave pairing: The peaks at (/> = (Af — l)/2 and (p = have 
the same sign, (b) d-wave pairing: The peaks at (p = {M — l)/2 and = have the opposite sign. 



defined on the lattice sites (s-wave), Eq. becomes 



M 



J = ^3c,i sin(%). 



(52) 



i=l 

If jc.i is equal for all one obtains a flux dependence similar to the Fraunhofer pattern: 

M^f ,.„(,,.)/(M + l)^,^"'W^^ + y,^'^), (53) 
J(0) ^ ^ ^ ' (M + l)sin(A;/M) ^ ' 

This formula reproduces the flux dependence of the supercurrent as obtained from the Bo- 
goliubov - de Gennes equations (shown in Fig. [T8|) . apart from slight deviations in the 
amplitude around the central peak at = (M — l)/2. It explains naturally the deviation 
from the $o/2 periodicity: it is an effect of discreteness, caused by the fact that the number 
of lattice sites in y-direction exceeds the number of plaquettes by one. 

/^From what has been explained for an s-wave junction, we construct a simple Ginzburg- 
Landau analogon for a rf-wave junction. In a d-wave superconductor, the order parameter 
is defined on the bonds between two neighboring lattice sites, and we therefore define the 
corresponding supercurrent as 

J= 5^j,,sin(A;(y, + l/2)). (54) 



For a constant jc.i we obtain 

J(0) 



M-\ 



Mk{y^ + l/2))/M 



sin (k) 



2=1 



Msin(A;/M)' 



(55) 
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FIG. 19: Current-flux relation calculated for a Josephson junction with Af = 11 and an inhomoge- 
neous impurity distribution. Top panel: The gray plaquettes in the profile of the junction have a 
scattering potential U = lOOt, while the white plaquettes have U = 2t, thus leaving two transparent 
channels through which almost the entire current flows. Blue: Bogoliubov - de Gennes equations 
and turquoise: Ginzburg-Landau model for (a) s-wave pairing and (b) d-wave pairing. 

which indeed reproduces the $o/2 periodic Fraunhofer pattern obtained from the Bogoliubov 
- de Gennes equations with nearest-neighbor pairing. The deviations in the amplitude are 
larger than for the s-wave junction telling that the d-wave junctions fulfill the Ginzburg- 
Landau conditions not as well as the s-wave junction. 

The Ginzburg-Landau formulae ( l52l) and ( IMI) are suitable also to calculate the supercur- 
rent flowing across junctions with an inhomogeneous impurity distribution. It is instructive 
to compare also in this case the supercurrent to results obtained from the Bogoliubov - de 
Gennes equations. Figure [19] shows such a comparison for a junction with M = 11 and 
current flowing only through the two "gaps" between the white plaquettes in the top panel 
of Fig. [191 In the microscopic model, this is achieved by setting strong repulsive poten- 
tials U = loot on the sites of the gray plaquettes and a small potential U = 2t on the 
white plaquettes. In the Ginzburg-Landau approach, we set jc,j = except for the two 
transparent channels. This system appears to be quite far from respecting the conditions 
for the validity of the Ginzburg-Landau equations. Nevertheless, for the s-wave junction, 
the results obtained from the Bogoliubov - de Gennes and Ginzburg-Landau equations are 
remarkably close. Even for the d-wave junction, the simple implementation of the Ginzburg- 
Landau equations reproduces the same features as the Bogoliubov - de Gennes equations, 
in particular it has maxima for similar values, but the amplitudes of the oscillations deviate 
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strongly. 

These considerations jointly lead to the conclusion that, even for small junctions where 
discreteness is pronounced, we do not find any indications that the Doppler shift has an effect 
on the current-fiux relation of Josephson junctions in the tunneling regime. The essential 
characteristics of the current-fiux relation, especially the position of the current maxima, 
agree quite well with the Ginzburg-Landau approach, where these effects are not included. 

2. Current-flux relation of transparent junctions 

A magnetic field threading a Josephson junction generates a supercurrent circulating 
around the junction, similar to a vortex in a type II superconductor, but with the complete 
fiux confined to the junction. If the junction is sufficiently transparent, the order parameter 
reacts to the current loop with a phase winding as in a fiux-threaded ring, with a winding 
number q that minimizes the total energy. The superconducting state in a transparent 
junction is therefore characterized similarly as a loop by the quantum number q related to a 
center-of-mass motion of the Cooper pairs and the supercurrent across the junction changes 
sign when the condensate reconstructs to another q. Remarkably, if the transparency is 
reduced, the discontinuities vanish smoothly, the current-fiux relation of the superconducting 
state with fixed q becomes periodic in 0, and in the tunneling regime, all states with different 
q become equivalent. This behavior is illustrated in Fig. [201 which shows E{(j)), J{(f)), and the 
spectrum for a uniform junction with nearest-neighbor pairing and U = 2t. The total energy 
consists of a series of parabolae, which correspond to different phase winding numbers. The 
kinks in E{(j)) and in the flux dependence of the spectrum are sharp for small values of 0, 
but the flnite repulsion on the junction smoothens the discontinuities in the supercurrent. 
Although the Doppler shift of the energy levels is not strongly pronounced in Fig. [201 the 
physical phenomena typical for multiply connected geometries govern the fleld dependence 
of the supercurrent across a Josephson junction, if its transparency is sufficiently high. 

V. CONCLUSIONS 

For unconventional nodal superconductors we established within a momentum-space for- 
mulation for superconducting loops that hc/e oscillations are present in the fiux dependence 
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FIG. 20: Characteristics of a transparent junction as obtained from the Bogohubov - de Gennes 
equations in a system with = 14 and M = 12 and a homogeneous impurity distribution with 
repulsive potential U = 2t. (a) The total energy, (b) the Josephson current, and (c) the energy 
spectrum of the system versus the external flux through the junction. 

of the ground state. The calculations in momentum space were restricted to rotationally 
symmetric systems like a cylinder, the energy spectrum of which depends sensitively on 
microscopic details. In Sec. 1111 Bl we have provided an analysis of the flux periodicity in a 
square frame with (i-wave pairing symmetry analogous to the cylinder geometry of Sec. [TTl 
with remarkably similar results. Nevertheless, the real-space calculations contributed to the 
understanding of the flux periodicity. We verified that the characteristic flux dependence of 
the (i-wave energy spectrum does not depend on the geometry or the absence of impurities. 

Within the real-space formulation, we constructed and analyzed more complex systems, in 
particular we investigated the periodicity of Josephson junctions. The idea that the Doppler 
shift drives energy levels through the Fermi energy in junctions between d-wave supercon- 
ductors, and thereby doubles the periodicity of the current-phase relation, seemed natural, 
but the physics turned out to be more subtle. Narrow junctions with only a few channels 
always display a period in the phase difference of An, even for s-wave superconductors, and 
the Doppler shift in tunnel junctions is too small to influence the current-phase relation. 
Only for transparent junctions does the Doppler shift become important; in this regime the 
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supercurrent across a Josephson junction behaves similar to the persistent supercurrent in 
a loop. These observations are also valid for the current-flux relation of fleld-threaded junc- 
tions. The microscopic theory excellently reproduced the results from the Ginzburg-Landau 
description of Josephson junctions in the tunneling regime, even for nanoscopically small 
systems with d-wave pairing. 

We thank R. Fresard, S. Graser, C. Schneider, and D. Vollhardt for insightful discussions and 
J. Mannhart for valuable contributions. This work was supported by the Deutsche Forschungsge- 
meinschaft through SFB 484 and TRR 80. 
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